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We present a comprehensive computational study of the phase diagram of the frustrated S = 1/2 Heisenberg 
antiferromagnet on the honeycomb lattice, with second-nearest ( J2) and third-neighbor ( J3) couplings. Using 
a combination of exact diagonalizations of the original spin model, of the Hamiltonian projected into the nearest 
neighbor short range valence bond basis, and of an effective quantum dimer model, as well as a self-consistent 
cluster mean-field theory, we determine the boundaries of several magnetically ordered phases in the region 
J2, J3 G [0, 1], and find a sizable magnetically disordered region in between. We characterize part of this mag- 
netically disordered phase as a plaquette valence bond crystal phase. At larger J2, we locate a sizable region 
in which staggered valence bond crystal correlations are found to be important, either due to genuine valence 
bond crystal ordering or as a consequence of magnetically ordered phases which break lattice rotational symme- 
try. Furthermore we find that a particular parameter-free Gutz wilier projected tight-binding wave function has 
remarkably accurate energies compared to finite- size extrapolated ED energies along the transition line from 
conventional Neel to plaquette VBC phases, a fact that points to possibly interesting critical behavior - such as a 
deconfined critical point - across this transition. We also comment on the relevance of this spin model to model 
the spin liquid region found in the half-filled Hubbard model on the honeycomb lattice. 

PACS numbers: 75.10.Kt, 75.10.Jm, 75.40.Mg 



Magnetic frustration is a very appealing route to weaken 
or destroy magnetic order, which can result in new phases of 
matter: these phases can usually be classified and named ac- 
cording to the broken symmetry (spin, lattice) if any, or they 
can belong to the spin-liquid zoo when no symmetry is bro- 
ken^ . The quest for a genuine gapped spin-liquid in a spin- 1/2 
model with SU(2) symmetry and an odd number of sites in 
the unit cell has started a long-time ago with the proposal by 
Anderson^ that the ground- state of the Heisenberg model on 
the triangular lattice could be viewed as a superposition of 
short-range valence bonds (VB), called a resonating valence 
bond (RVB) state. For the specific example of the triangular 
lattice it turned out later however that a magnetically ordered 
state is realized^. Up to now, there is still no firmly estab- 
lished spin-liquid ground-state with the aforementioned prop- 
erties in a reasonably realistic SU(2) spin model, although 
there are potential candidates, such as the triangular lattice 
with ring exchange interactions^ or the Heisenberg model on 
the kagome lattice^ . On the other hand, if one considers lat- 
tices with an even number of sites per unit cell, then Hasting 's 
generalization^ of the Lieb-Schultz-Mattis theorem^ does not 
apply, and it is possible in principle to stabilize a magnetically 
disordered ground- state that does not break any symmetry and 
has only trivial topological properties. One can think for in- 
stance of a Heisenberg model on a square bilayer lattice with 
strong interlayer exchange. The honeycomb lattice is peculiar 
in this respect because no simple lattice- symmetry preserving 
deformation is known which would lead to a gapped magnetic 
stated 
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FIG. 1 . (Color online) (a) Honeycomb lattice with the different spin 
exchange interactions considered in this paper; (b) corresponding 
Brillouin zone with relevant k points. 



In recent years a promising new direction in the search for 
spin liquids has opened up, focusing on the behavior of insu- 
lating phases upon approaching the Mott insulator-metal tran- 
sition. In the half-filled triangular lattice Hubbard model a 
picture with a spin bose metal spin liquid phase sandwiched 
between the metallic phase at small U/t and the magnetically 
ordered Neel phase at large U/t has emerged^ It has been 
recognized that this spin liquid phase can be understood in 
terms of a pure spin model, where the rising charge fluctua- 
tions are cast into an increasingly complex spin Hamiltonian 
beyond the Heisenberg model^^'^"^. A second striking exam- 
ple of a spin liquid located between a magnetically ordered 
phase and a (semi)metal has recently been uncovered in the 
half-filled Hubbard model on the honeycomb lattice^^. Such 
spin liquid phase is reported to have a small spin gap and no 
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FIG. 2. (Color online) Phase diagram of the frustrated S = 1/2 Heisenberg model honeycomb lattice in the region J2, Js ^ [0, 1], based on a 
combination of exact diagonalization results discussed in the main text. The 5 regions identified here correspond to: (I) a Neel ordered phase 
with staggered magnetization, (II) a collinear magnetically ordered phase, (III) One or several phases corresponding to short or long range 
ordered non-collinear magnetic order, (IV) A different collinear magnetically ordered (or disordered) phase corresponding to phase (IV) in 
Ref. 9 and (V) a magnetically disordered phase forming a plaquette valence bond crystal. The five phases are sketched in the panels around the 
phase diagram. Note that the phases highlighted in grey (III), (IV) show substantial finite size effects and are therefore difficult to characterize 
precisely. 



appreciable correlations of any kind. 

This exciting finding leads us to the natural question of 
whether this spin liquid phase on the honeycomb lattice can 
also be described within a pure S = 1/2 spin model, de- 
spite the vicinity of the insulator to semimetal transition. A 
high order derivation of the corresponding spin model is in 
progress but the typical value of the expansion parame- 
ter t//7 ~ 0.25 relevant for the spin liquid phase renders 
this task more challenging in comparison to the triangular lat- 
tice, where a typical value for the spin liquid regime is about 
t/U ~ 0. 1 1 . In the absence of an accurate prediction for a rel- 
evant spin model, we start by exploring the effect of the next- 
to-leading order correction to the nearest neighbor Heisenberg 
model, which is a second neighbor Heisenberg coupling J2 
arising at fourth order in t/U . We thus consider in the fol- 
lowing a frustrated 5 = 1/2 Heisenberg Hamiltonian on the 
honeycomb lattice, where we also include a third neighbor 
coupling J3 for completeness. 

The honeycomb (hexagonal) lattice Hamiltonian [see 
Fig. 1(a)] reads: 

(m) {{i.j)) {{{hj))) 



In this paper we focus solely on antiferromagnetic interac- 
tions Ja > 0, set Ji = 1 and restrict ourselves to the win- 
dow J2, J3 G [0, 1]. Aspects of this frustrated model have 
been explored previously in the literature, based on spin wave 
theory^' a nonlinear sigma model treatment^^, Sch winger 
boson approaches^ ^'^^ and exact diagonalizations^'^^ . Note 
also that a similar frustrated model gives rise to a rich phase 
diagram on the square lattice^"^'^^ . 

In this work we thoroughly explore the phase diagram in 
the considered window, based on a combination of exact di- 
agonalizations (EDs) of the spin model (up to 42 spins), EDs 
in the nearest-neighbor valence bond (NNVB) subspace (up to 
96 spins), EDs of an effective quantum dimer model (QDM) 
(corresponding to up to 126 spins), complemented by a self- 
consistent cluster mean field theory (SCMFT) and the study 
of a fully projected Gutz wilier wave function of the half-filled 
honeycomb tight binding "Dirac sea". 

The key finding of our work is the presence of a sizable 
magnetically disordered region adjacent to the well studied 
Neel phase of the unfrustrated honeycomb Heisenberg model. 
We identify a large part of this region as a plaquette valence 
bond crystal (VBC). Interestingly we find evidence (within 
the ED realm) for a possibly continuous phase transition be- 
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tween the Neel phase and a plaquette valence bond crystal. 
In addition the energy and some of the key correlations of the 
frustrated spin model in the transition region are well captured 
by a simple Gutz wilier projected (GP) "Dirac sea" wave func- 
tion. These findings raise the possibility of a continuous quan- 
tum phase transition beyond the Ginzburg Landau paradigm in 
this honeycomb lattice spin model. 

The outline of the paper is as follows: we start by giving 
a quick overview of the phase diagram in section I. Then the 
magnetically ordered phases are located using a SCMFT in 
section II. Next we study the spin model using EDs in sec- 
tion III, followed by EDs in the NNVB subspace and EDs 
of a QDM which are presented in section IV. We close with 
a discussion and conclusion section V. In the appendices we 
discuss the properties of the Gutzwiller projected "Dirac sea" 
(App. A), present the derivation of an effective quantum dimer 
model from the Hamiltonian projected into the NNVB sub- 
space (App. B), compare the energies and the finite size be- 
havior of the NNVB versus the QDM approach (App. C) and 
derive the expected correlation functions in model valence 
bond crystal states (App. D). 



pled triangular lattices in the large J2 limit; (IV) a collinear 
magnetically ordered phase (corresponding to phase (IV) in 
Ref . 9) or a staggered dimer phase (also called lattice nematic 
in Ref. 18). And finally (V) a magnetically disordered phase 
forming a plaquette valence bond crystal. 

Analyzing the magnitude of the fidelity dips it seems likely 
that the transition from phase (I) to (V) is continuous (corre- 
sponding to a faint feature in the fidelity), while the transitions 
from (I) to (II) and (V) to (II), (III) and (IV) seem to be of first 
order because of strong avoided level crossings observed on 
the clusters considered. The topology of the phase diagram 
and the nature of the phases in the regions (III) and (IV) dis- 
play strong finite size effects and require further investigations 
beyond the scope of this work. We note that phases (II), (III) 
and (IV) exhibit long range staggered dimer order also in the 
case of long-range magnetic order, because in the magneti- 
cally ordered phases one of the nearest neighbor (NN) bond 
energies is different from the other two. 

We now proceed to a self consistent cluster mean field treat- 
ment which is well suited to detect various magnetically or- 
dered phases. 



I. OVERVIEW OF THE PHASE DIAGRAM 

We start by summarizing the main result of this paper, the 
phase diagram of the frustrated 5 = 1/2 Heisenberg Hamil- 
tonian Eq. (1) in the considered parameter window displayed 
in Fig. 2. The phase diagram emerges from a combination of 
different information extracted from ED of the spin model: 

(i) For a first analysis without further input we have inves- 
tigated the structure of the fidelity /, i.e. the overlap between 
ground-states (GS) obtained for different parameters'^: 

f{J2,MJ2,Jk) = \{GS{J2,J3)\GS{J^,4))\, (2) 

for consecutive points along the two directions of the (72,^3) 
plane using a grid spacing of 0.05. Local minima of / in the 
directions of J2 and J3 of / are indicated by star symbols for 
both the N = 24 and N = 32 samples. These dips already 
give a first impression of some phase boundaries in the phase 
diagram. 

(ii) In addition we highlight the quantum numbers of the 
lowest excited state and - if not a triplet already - the quantum 
numbers of the lowest triplet, both for N = 24. The basic idea 
is that for sufficiently large systems the quantum numbers of 
the low energy spectrum are characteristic of the respective 
phases, and can thus be used to chart a phase diagram if used 
with care. For a detailed discussion of the expected low ly- 
ing energy levels in the different phases we refer to subsec- 
tion III A. 

Based on these and further results to be discussed later, the 
following phases are identified: (I) a Neel ordered phase with 
a finite staggered magnetization, located around the unfrus- 
trated point J2=J3=0; (II) a magnetically ordered collinear 
phase corresponding to the classical phase (II) in Ref. 9 arising 
at combined large J2 and J3 ; (III) one or several phases cor- 
responding to short-range or long-range non-collinear mag- 
netic order, resulting from the Ji , J3 coupling of two decou- 



II. SELF CONSISTENT CLUSTER MEAN-FIELD 
THEORY 

The very same frustration accounting for the rich physics 
exhibited by the here considered model [Eq. (1)] also adds 
enormous complexity to the task of determining its properties. 
In this context, approximate approaches can be valuable and 
employed in obtaining "draft phase diagrams" that may guide 
subsequent application of more accurate techniques. The so- 
called self-consistent cluster mean-field theory (SCMFT)^^'^^ 
is a tool particularly well suited to this task. In a nutshell, 
SCMFT consists in diagonalizing the Hamiltonian under in- 
vestigation on small clusters that, besides including actual 
in-cluster couplings, so that quantum fluctuations at the lo- 
cal level are partially taken into account, are also coupled to 
mean fields that are to be determined self-consistently. This 
technique has been shown to considerably improve upon more 
conventional mean-field approaches for the case of hard-core 
bosons on the triangular lattice^^ and, more recently, to yield 
results that compare well with the ones from more sophis- 
ticated techniques when applied to an effective model for a 
frustrated antiferromagnet.^^ 

In applying SCMFT, we consider the clusters comprising 
N = 6 and 8 sites depicted in the insets of Fig. 3. We split the 
Hamiltonian Eq. (1) according to 

?i=?iin+?iMF. (3) 

Hin accounts for in-cluster couplings, 
7iin = ^S,-S, + J2 ^ S,-S,+J3 Yl S.-S,,(4) 

and is treated in an exact way. (i, j), {{hj)) and (((i, j))) re- 
spectively denote nearest-, second-nearest- and third-nearest- 
neighbor in-cluster sites (open circles in Fig. 3, where in- 
cluster NN bonds are represented by thick lines). Couplings 
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FIG. 3. (Color online) Phase diagram obtained from SCMFT as applied to the = 8 (a) and = 6 (b) clusters. In the region labeled as 
"disordered" no S*?/ (2) -broken magnetic phases are obtained from the SCMFT procedure for the considered cluster. Insets: clusters employed 
in our SCMFT calculations. Thick lines connecting open circles represent in-cluster couplings and dashed lines coupling to mean fields (only 
Ji interactions are depicted). 



to the mean fields are included in Hmf that reads 



n 



MF 



= ^S,-(S,) + J2^S,-(S,)- 



-Js 

[[[iJ]]] 



(5) 

Here, denotes a spin-operator attached to an in-cluster site i 
and is coupled to the mean-field given by the expectation value 
(Sj) at the "across-the-boundary site" j [light-filled circles 
in Fig. 3], for nearest ([z, j]; dashed lines in Fig. 3), second- 
nearest ([[i, j]]) and third-nearest ([[[^, j]]]) neighbors. That 
is, one may see SCMFT as an exact diagonalization on a fi- 
nite cluster with periodic boundary conditions (PBC), where 
"across the boundary" interactions are replaced by couplings 
to mean fields that are determined in a self-consistent manner. 
One starts from a randomly chosen wave function and com- 
putes the mean fields (Sj) at every site j, that are then used 
in setting Hmf- The cluster Hamiltonian Eq. (3) is then di- 
agonalized and the so-obtained ground- state wave function is 
used in re-setting Hmf; computation proceeds until all (Sj) 
are converged and the existence of 6'/7( 2) -broken magnetic 
phases is signaled by non- vanishing mean fields, (Sj) ^ 0. 

In Fig. 3 we present the resulting SCMFT phase diagrams 
obtained for (a) an = 8 cluster and (b) an A/" = 6 clus- 
ter^^. For the A/" = 8 cluster we first note the presence of two 
collinear magnetically ordered phases, labeled (I) and (II) in 
Fig. 2. These phases are also present in the classical version 
of the model, occupying roughly the same portion of the plane 
{J2j Js)^'^^ • Furthermore, the phase labelled as IV in Fig. 2 — 
also observed in the classical case but only for J3 < 0,^'^^ — 
occupies part of the region shown to support a spiral phase 
in Ref. 17. This might be an interesting effect, where the 
collinear phase IV is stabilized for some J3 > (i.e. be- 



yond the classical domain of stability) by quantum fluctua- 
tions. Note that a magnetically ordered phase of this type is 
also compatible with the pronounced staggered dimer pattern 
reported in previous ED studies^'^^ for J2 > 0.4, and the lat- 
tice nematic point of view^^. 

In order to study finite size effects we apply SCMFT to an 
N = 6 site cluster and present its phase diagram in Fig. 3(b). 
First, we remark that this cluster [depicted in the inset of 
Fig. 3(b)] is not compatible with both phases II and IV and that 
no solutions with (Sj) 7^ are encountered in some parts of 
the region stabilizing these orderings for the A^ = 8 site clus- 
ter [Fig. 3(a)]. Furthermore, the size of the region supporting 
Neel order is somewhat reduced^ ^ in comparison with what is 
observed in Fig. 3(a): we suspect that this may be explained 
by the fact that the Kekule-like state with resonating valence 
bonds is particularly stable on the hexagon- shaped A^ = 6 
cluster, with the consequence that the "disordered" region is 
overestimated. More interestingly, and in contrast with what 
happens for the A^ = 8 cluster, a spiral state (phase III) is sta- 
bilized for large J2 . Such a state is adiabatically connected to 
the ground-state for J2/J1 ^ 1, where Eq. (1) decouples into 
two triangular lattices, each of which exhibits 120° magnetic 
order. 

Finally, we address the possible occurrence of non- 
magnetic phases for Eq. (1). We notice the existence of an 
extended region in (72,^3) where vanishing mean-field solu- 
tions, (Sj) = 0, are obtained for both clusters considered in 
Fig. 3. Intersecting the magnetically disordered phases from 
both phase diagrams (N = 6 and A^ = 8) we obtain a putative 
nonmagnetic region which roughly corresponds to the extent 
of phase V in the ED based phase diagram shown in Fig. 2. 
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In the following we address the nature of this magnetically 
disordered region in more detail and determine some of the 
phase boundaries with higher accuracy based on finite size 
extrapolated ED simulations. 

III. EXACT DIAGONALIZATION 

We now explore the phase diagram based on large scale 
ED in the basis of finite honeycomb lattice samples with 
N = 24,26,28,30,329'23,34,36,38^^ and 42 sites. The 
clusters with N = 24,30,36,42 sites feature the two K 
points in their Brillouin zone [c.f. Fig. 1(b)], while = 
24, 28, 32, 36 contain one or several M points. The clusters 

= 24, 26, 32, 38, 42 exhibit six-fold rotational symmetry. 

We first study the low-lying energy spectrum in the full pa- 
rameter region in subsection III A in order to provide more 
information on the phase diagram shown in Fig. 2. Then we 
address the stability of the Neel phase (I) by calculating mag- 
netic structure factors and energy scalings in subsection IIIB, 
and close this section with a discussion of the nature of the 
dimer-dimer correlations in subsection IIIC, supporting the 
presence of an extended plaque tte valence bond crystal phase. 

A. Nature of the lowest excitation 

In compiling the phase diagram shown in Fig. 2, one set 
of information was gathered from the quantum numbers of 
the low-lying excitations. The idea is that symmetry broken 
phases must exhibit a specific set of low lying energy levels, 
which will allow the spontaneous symmetry breaking in the 
thermodynamic limit. In the case of SU {2) symmetry break- 
ing states, the appropriate structure is called "tower of states" 
(TOS) and has been successfully used to identify magnetically 
ordered phases^, as well as spin nematic phases^"^. The finite 
size behavior of energy gaps is as follows: the levels belong- 
ing to the symmetry breaking tower states scale as 1/7V with 
system size, while the spin-wave modes scale as 1 /L ^'^^ . For 
sufficiently large system sizes one should therefore detect only 
states belonging to the TOS manifold in the lowest part of the 
energy spectrum. In the case of discrete symmetry breaking 
- such as for a VBC - a finite number of levels is expected to 
collapse rapidly (exponentially beyond a certain correlation 
length) onto the ground state. In each case, the quantum num- 
bers of the collapsing levels are determined by the nature of 
the order parameter, i.e. the broken symmetries, and generally 
they will be different for distinct phases (but not always). In 
the following we summarize the expected quantum numbers 
of low energy levels of several candidate phases (some of the 
results were presented earlier in Ref. 9). Note that the quoted 
quantum numbers are given as appropriate for the = 24 
sample, and the C^y point group is located at the center of a 
hexagon. 

1 . The Neel ordered phase (I) has a simple TOS structure 
with one level per total spin: all even spin sectors be- 
long to the F Al representation, while the odd ones ap- 
pear in F B2. 



2. The magnetically ordered phase (II) has three levels per 
spin sector: the levels in the even spin sectors are found 
in F Al and the two-dimensional representation F E2. 
The levels in the odd spin sectors belong to the threefold 
degenerate M momentum, with even (odd) parity for 
reflections along (perpendicular) to the F — M axis. 

3. The magnetically ordered phase (IV) has three levels 
per spin sector: the levels in the even spin sectors are 
found in F Al and the two-dimensional representation 
F E2. The levels in the odd spin sectors belong to the 
threefold degenerate M momentum, with even parity 
for both reflections along and perpendicular to the F — 
M axis. 

4. A columnar (Read-Sachdev^^) [cf. Fig. 4 (a)] or pla- 
que tte VBC [cf. Fig. 4 (b)] has three collapsing singlet 
levels: one at F Al and a two-dimensional K Al repre- 
sentation. Note that these two VBCs can not be distin- 
guished based on energy level quantum numbers alone. 

5. A staggered VBC [cf. Fig. 4 (c)] has three collapsing 
levels: F Al and F E2 (2-dim representation). 

It is interesting to note that the level crossing of the excited 
states quantum numbers shown in Fig. 2 match quite well the 
dips in the fidelity (which is a ground state observable). The 
quantum numbers carry however more information and allow 
to label the quantum phases roughly, before studying them in 
more detail using correlation functions, as we will do in the 
following. 

B. Stability of the Neel phase 

The non-frustrated model (J2=-^3=0) is known to pos- 
sess antiferromagnetic (AF) long-range order. This has been 
shown by several techniques including linear spin-wave the- 
ory^^, a coupled cluster method^^, ED^^'^^, series expansions 
around the Ising limit^^, tensor network studies^^"'^^ varia- 
tional Monte-Carlo"^^ and quantum Monte-Carlo (QMC) sim- 
ulations"^^""^^ . In particular, the staggered moment is m^o = 
0.2677(6)'^'^, a value that is significantly reduced by quantum 
fluctuations compared to the classical value of 1/2. 

In Fig. 5, we plot ED data for the finite- size magnetic order 
parameter squared"^^: 

for various clusters sizes N and J2 values (we set J3 = 
for the moment). Standard finite-size scaling predicts leading 
1/L = 1/ \/N corrections^^'^^, which we find to be quite well 
satisfied even for small clusters for the unfrustrated case. The 
infinite system size estimate including all system sizes shown 
(N = 24, 26, 28, 30, 32, 34, 38) is m'^{oo) = 0.0684, (rrioo = 
0.262)^^. Our best agreement with QMC is found based on 
the samples with TV = 24, 28, 32 sites only, yielding an esti- 
mate of m^(oo) = 0.0728, corresponding to moo = 0.270. 
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FIG. 4. (Color online) Pictorial representation of the three valence bond crystal candidate states discussed in this work. The columnar VBC 
(a) is also called Read-Sachdev^^ state in the literature, while the staggered dimer VBC is also known as "lattice nematic"^^. 
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FIG. 5. (Color online) Squared staggered moment vs 1/L at 
Js = for several J2 values obtained by ED on the clusters N = 
24, 26, 28, 30, 32, 34, 38 and the corresponding extrapolations to the 
thermodynamic limit. In the non-frustrated case J2 = 0, we obtain 
a good agreement with the QMC value when using the N = 24, 28 
and 32 samples alone, (see text). Inset: extrapolated value of the 
staggered moment as a function of J2/J1, vanishing between 
^2/^1^ 0.17 — 0.22, depending on the extrapolation. 



The discrepancy between the different ED extrapolations re- 
sults in a ~ 5% uncertainty on the magnitude of the magnetic 
moment. 

When J2 is switched on we notice that the finite size data 
starts to deviate systematically from a straight line in We 
observe that systems which contain an M point in the Bril- 
louin zone (N = 24, 28, 32), behave consistently with respect 
to each other - studying e.g. the derivative dm^/dJ2 - com- 
pared to the other system sizes"^^. We therefore choose to base 
one of the extrapolations (solid circles in Fig. 5) on this class 
of samples. The second estimate is obtained by using all the 
shown system sizes (hatched circles in Fig. 5). Now, as J2 is 
increased starting from zero, the extrapolated staggered mo- 
ment m^(oo) decreases quite rapidly, roughly linear with in- 
creasing J2 (inset of Fig. 5), and vanishes continuously around 
J2 = 0.17 ~ 0.22, based on the two extrapolations. Despite 



some uncertainty, this constitutes a critical value which is 
larger than the classical estimate 1/6^^, the linear spin wave^ 
and non-linear sigma model^^ results of 0.1 ~ 0.12, and sub- 
stantially larger than a recent variational Monte Carlo (VMC) 
estimate of 0.08^^. Our estimate is however in agreement with 
a Sch winger-boson mean field treatment which reported a crit- 
ical J2 of about 0.2^^. A possible physical explanation of this 
shift of the transition to larger values of J2 is that in some 
cases quantum fluctuations prefer collinear over spiral states, 
such as e.g. in the Ji — J3 model on the square lattice^'^'^^'^^ 

We have also determined the ordered moment along a sec- 
ond J2 cut at constant J3 = 0.3 (data not shown). In this 
case the cluster size and shape dependency is even more pro- 
nounced and makes an accurate determination of the critical 
J2 value rather difficult. Similar extrapolations based on ei- 
ther all samples or only a subset of the samples gives a tran- 
sition point somewhere between J2 ~ 0.27 — 0.33, although 
the actual uncertainty is probably larger. 

In order to corroborate the location of the disappearance of 
Neel order, we study the energy per site along the same two 
constant J3 lines, one located at J3 = and the other one 
at J3 = 0.3. In a Neel ordered phase the leading finite size 
corrections to the energy per site are expected to scale as 



E/N = eN = ec 



a c 



(7) 



i.e. with a leading correction^^'^^, and the coefficient of 

this term is proportional to the spin wave velocity c. In Fig. 6 
(J3 = on the left panel and J3 = 0.3 on the right panel) 
we display the energy per site of samples of up to 42 spins 
in the relevant J2 range together with the resulting ^ 00 
estimate Coq. We extrapolate the energy according to Eq. (7) 
up to J2 = 0.3 for J3 = 0, and up to J2 = 0.4 for J3 = 0.3. 
As can be seen in Fig. 7, the prefactor of the cor- 
rection term is reduced upon approaching the transition re- 
gion, but seems to stay constant at the transition, in analogy 
to the frustrated square lattice antiferromagnet"^^'^^. Note that 
the shaded regions denote the approximate locations of the 
transitions based on the extrapolation of the ordered moment, 
and the minimum of the velocity agrees reasonably well with 
those estimates. 

Returning to the extrapolated energies, we note that for the 
unfrustrated J2, J3 = case in Fig. 6 (left panel) the extrapo- 



7 




FIG. 6. (Color online) Left Panel ED: Ground state energy of samples up to = 42 at J3 = 0, together with the resulting N ^ 00 estimate. 
For comparison the QMC result at J2 = ^^'"^^ and the energy of the Gutzwiller projected Dirac sea (see App. A) are displayed. The light red 
shaded area denotes the approximate location of the disappearance of Neel order, while the dark blue shaded region denotes the approximate 
location of the first order transition to yet another phase (likely phase IV). Right Panel ED: Ground state energy of samples up to = 38 at 
J3 = 0.3 Ji, together with the resulting N ^ 00 estimate. For comparison our ALPS QMC result at J2 = and the energy of the Gutzwiller 
projected Dirac sea (see App. A) are displayed. The light red shaded area denotes the approximate location of the disappearance of Neel order, 
while the dark blue shaded region denotes the approximate location of the first order transition to yet another phase (likely phase III). Note 
for both panels the good agreement between the extrapolated ED and QMC results at J2 = on the one hand, and between the Gutzwiller 
projected Dirac sea wave function and the extrapolated ED data at the magnetic to non-magnetic transition on the other hand. 
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FIG. 7. (Color online) J2 dependence of the correction term 

to the energy per site. This quantity is proportional to the spin wave 
velocity c in the Neel ordered phase. At the transition its value is 
expected to be finite. 



lated energy per site e 00 is in very good agreement with pub- 
lished QMC results^^'^^, and we find similarly good agree- 
ment for J2 = 0, J3 = 0.3 in Fig. 6 (right panel), where we 
performed ALPS SSE simulations^^'^"^ to obtain an accurate 
estimate for the energy. For both J3 values the energy then 
first rises almost linearly with increasing J2 , as expected for 
this particular Neel phase (note that the derivative de/dJ2 is 
proportional to (S^ • Sj) on the J2 bonds as a consequence of 
the Hellmann-Feynman theorem). The energy curves flatten 



at larger J2 and exhibit a maximum around J2 ~ 0.35 — 0.4 
for J3 = and J2 ~ 0.5 — 0.6 for J3 = 0.3. A comparison 
with the fidelity data shown in Fig. 2, suggests that the maxi- 
mum of the energy approximately coincides with the avoided 
level crossing to a different phase. 

Inspired by the success of a simple Gutzwiller projected 
half-filled tight-binding wave function on the triangular lat- 
tice in describing the spin liquid regime on the insulating side 
of the Mott transition we have analyzed a related wave 
function on the honeycomb lattice: the Gutzwiller projected 
half-filled honeycomb tight binding wave function (termed 
the Gutzwiller projected "Dirac sea" in the following). This 
wave function is discussed in some detail in App. A. It is a 
parameter-free variational wave function, and its energy for 
the Hamiltonian considered here is given in Eq. (Al). This 
energy is plotted using a dashed green line in Fig. 6. Quite 
remarkably the energy is very close to the finite size extrapo- 
lated ED energies precisely in the region where the Neel or- 
der is about to vanish (light red shaded uncertainty regions). 
As discussed in App. A, the Gutzwiller projected Dirac sea 
wave function has algebraically decaying spin- spin correla- 
tions with a sign structure that is identical to the one displayed 
by the Neel state. With its algebraically decaying correlation 
this wave function could in principle describe qualitatively a 
putative continuous quantum phase transition from the Neel 
ordered phase to a quantum paramagnet. Inspecting the na- 
ture of the dimer-dimer correlations in the Gutzwiller pro- 
jected Dirac sea, the signs of the dimer-dimer correlations are 
identical to the ones in the columnar (Read-Sachdev) or pla- 
quette VBC (disregarding one particular distance highlighted 
in Fig. 20). Given the surprisingly accurate energy of this 
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FIG. 8. (Color online) Four-spin correlations [Eq. (8)] on an = 24 site cluster for several different J2 values and constant J3 = 0.3: (a) 
(0.0, 0.3) in the Neel phase I, (b) (0.3, 0.3) in region of the phase transition from phase I to V, (c) (0.5, 0.3) in the plaquette phase V, and 
(d) (0.7, 0.3) in the lattice nematic {staggered dimer VBC) phase III or IV. The reference bond is indicated by the thick-black line. Negative 
(positive) correlations are represented by red (blue) bonds. 



wave function at the transition, a plausible scenario is the pres- 
ence of a continuous Neel to columnar/plaquette VBC quan- 
tum phase transition in this frustrated honeycomb antiferro- 
magnet. We will discuss this scenario and other possibilities 
in more detail later on. Let us now consider the dimer-dimer 
correlations in the magnetically disordered phase, in order to 
verify whether there is indeed a VBC phase present. 

C. Dimer correlations 

Real space correlations. In this section we study dimer- 
dimer correlations using ED. Our aim is to highlight the struc- 
ture of the correlations along a J2-cut at constant J3 = 0.3 for 
the A/" = 24 sample^^. We measure the following four-spin 
correlation function: 

C^jki = 4 (((S, • S,) (Sfc • Si)) - ((S, • S,))') , (8) 

where i^j and k^l are nearest-neighbor bonds on the honey- 
comb lattice. In Eig. 8 correlation function results for four dif- 



ferent values of J2 are shown. The panel (a) shows the dimer- 
dimer correlations deep in the Neel phase at J2 = 0, J3 = 0.3, 
where we expect the correlations to decay rapidly, but with a 
power-law, due to the coupling to the multi spin wave con- 
tinuum. In panel (b), at J2 = 0.3, J3 = 0.3, we sit ap- 
proximately at the Neel to paramagnet transition, and some 
of the more distant bonds have changed sign compared to the 
Neel phase. Note that this correlation pattern matches qualita- 
tively the one of the Gutzwiller projected Dirac sea discussed 
in App. A (as well as the Z2 -liquid discussed in Ref. 50), and 
surprisingly also the one reported for the spin liquid regime 
of the half-filled Hubbard model in Ref. 15. Panel (c) at 
J2 = 0.5, J3 = 0.3 shows pronounced and long-ranged cor- 
relations, which at first sight seem to be compatible with ei- 
ther a columnar (Read-Sachdev) or plaquette VBC according 
to App. D. A more quantitative inspection reveals however 
that the largest distance positive/negative correlations are re- 
spectively close to 0.18 and -0.0935, which is in favor of a 
(d-wave) plaquette phase. We revisit the question of colum- 
nar versus plaquette order again in the context of the effective 
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FIG. 9. (Color online) Left panel: plaquette/columnar VBC structure factor 'S'p^^^ /coi./^b for N — 24, obtained using ED in the basis, 
as a function of ( J2, J3). The radius of the circles is proportional to *S'piaq./coi./^b. Numbers correspond to lO^S'pi^.q /Coi./^b- Right panel: 
staggered VBC structure factor ^SstagV^b for N — 24, obtained using ED in the basis, as a function of ( J2, J3). The radius of the circles 
is proportional to *S'stag./^b- Numbers correspond to lO^S'stagV^b- Note that the strong plaquette/columnar signal is found within phase 
(V) of Fig. 2, while the strong staggered signal is associated to the phases (II),(III) and (IV). 



quantum dimer model, and corroborate the present finding of 
a plaquette phase. Finally, panel (d) at J2 = 0.7, J3 = 0.3 
shows very strong correlations reminiscent of a staggered 
dimer phase. We stress again that this finding alone does not 
discriminate between a spin gapped valence bond crystal or a 
magnetically ordered phase of type (II), (III) or (IV). While 
a staggered valence bond crystal neighboring the plaquette 
phase is likely (according to Refs. 9, 18, and 23), at larger 
J2 the staggered signal in the dimer-dimer correlations could 
persist despite the appearance of magnetic order. 

VBC structure factors. It is instructive to grasp correlations 
using integrated quantities such as dimer structure factors. As 
displayed in Fig. 24 of App. D, two different dimer correlation 
patterns emerge for plaquette/columnar states or a staggered 
state. In order to detect them, we define the dimer structure 
factors as: 

>Sa = ^^a(^,OC^.fc^, (9) 

with a = Plaq./Col. ox a = Stag., where ^^(/c, /) = +1 if 
(/c, /) are NN sites such that Cijki > for "pure" a states and 
s^VBc(^,0 = ~2 otherwise (strong correlations the closest 
to the reference bond are not included for the related quantity 
5** ; see Ref. 24). It is important to stress here that the are 
order parameters detecting lattice symmetry breaking, which 
do not distinguish themselves between genuine VBC ordering 
or a lattice symmetry breaking magnetic state. 

A full scan of dimer structure factors associated to either 
plaquette/columnar or staggered valence bond crystal order 
is shown in Fig. 9. Consistently with real space dimer cor- 
relation analysis, two phases come up with strong plaque- 
tte/columnar (left panel) or staggered (right panel) signal. 




FIG. 10. (Color online) version of the dimer correlations obtained 
using ED in the basis for AT = 42 at J2 = 0.3, J3 = 0. 



The staggered signal is especially strong in the vicinity of 
the J3 = line and close to the avoided level crossing. In 
contrast, the columnar/plaquette signal is strongest around 
J 2 ~ 0.6, J3 ~ 0.4, and decreases upon approaching the 
J3 = line. In order to address the behavior at J3 = we 
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have also calculated the dimer correlations [cf. Eqn. (A2)] 
for J 2 = 0.3, J3 = on the N = A2 sample, which would be 
compatible with a columnar/plaquette VBC. The correspond- 
ing plot shown in Fig. 10 exhibits a correlation pattern remi- 
niscent of the one expected for columnar/plaquette states, but 
the correlations are not particularly strong, and also exhibit a 
few defects in the form of bond which show inverted corre- 
lations compared to the columnar/plaquette expectations. We 
are thus currently unable to discriminate whether this picture 
corresponds to a columnar/plaquette VBC with a small order 
parameter or a genuine spin liquid, and more work is needed 
to clarify the behavior at J3 ~ 0. 



IV. EXACT DIAGONALIZATIONS IN THE 
VALENCE-BOND BASIS AND QUANTUM DIMER MODELS 



A. Diagonalization in the Nearest-Neighbor Valence-Bond 
Basis 



In this subsection we present results obtained from exact 
diagonalizations in the (variational) basis given by the set of 
nearest-neighbor valence-bond states^"^'^^. Recently Mosadeq 
et alP presented an analysis using the same technique but 
limited to J3 = and small system sizes {N = 54). Here we 
consider the more general case of finite J3 G [0, 1] and inves- 
tigate considerably larger clusters (N = 72 for correlations 
and N = 96 for energies), allowing us to perform systematic 
finite size extrapolations. 



1. The Method 

When frustration is dominant and destabilizes magnetic 
phases it is possible to explicitly take into account that states 
with non-zero total spin are unimportant in accounting for the 
low-energy physics and to describe the system solely in terms 
of the 5 = subspace. This subspace can be spanned by the 
set of arbitrarily ranged valence-bond (VB) states^^'^^, which 
forms an over-complete basis and is thus difficult to manipu- 
late, especially in numerical studies. A natural way of circum- 
venting this difficulty is to impose a cutoff on the maximum 
range of the VBs to be considered; in particular, it is possi- 
ble to devise an approach where only nearest-neighbor VB 
(NNVB) states are taken into account. ^"^'^^ While the restric- 
tion to NNVB states is obviously a variational approximation, 
it offers the key numerical advantage of a significant reduction 
of the Hilbert space and has been shown to yield sound results 
for a number of strongly frustrated models, whose low-energy 
physics is dominated by short-range spin- singlets. ^"^'^^'^^ 

We briefly recall how the method can be applied and refer 
to Refs. 24 and 56 for details. We follow a heuristic argu- 
ment and try to formulate the eigenvalue problem in the re- 
stricted NNVB subspace simply as: X]^Q^iH|v^i) = 
E^-ai\(pi). However, since the set is not invariant 
under the application of H, this relation cannot hold in the 
particular singlet subspace but can explicitly be enforced in 
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FIG. 1 1 . (Color online) Comparison between the ground- state energy 
obtained from diagonalization in the NNVB and in the basis, as a 
function of J2, J3 G [0, 1]. The radius of the circles is proportional 



to the percent relative error, (Eq^^^ 
obtained from a cluster comprising N - 



- ^0 )/^o X 100. Data 
24 sites. 



the restricted NNVB subspace by considering 

^ai{(pj\n\(pi) = E^ai{(pj\(pi) , (10) 

for all \(pj) G This last equation is nothing but a gen- 

eralized eigenvalue problem (GEP) for the two matrices with 
elements given by Hij = {(pj\ii\(pi) and Oij = {(pj\(pi), 
the latter explicitly denoting the non-orthogonality of NNVB 
states. However, it is crucial here that despite of their non- 
orthogonality, the NNVB states are linearly independent on 
most relevant lattices^^'^\ and particularly on the honeycomb 
lattice with periodic boundary conditions considered here^^ 
GEPs are computationally more demanding than conventional 
eigenvalue problems, especially in the present case where both 
Hij and Oij are dense matrices. In spite of this, since for a 
given system size the dimension of the NNVB subspace is 
much smaller than that of the total = subspace, the 
method discussed here allows us to treat considerably larger 
clusters, and thus to perform more extended finite size extrap- 
olations, than possible within conventional ED. We also re- 
mark that, as for standard ED, it is possible to take advantage 
of lattice symmetries, so that the size of the matrices to be con- 
sidered is further reduced and, even more importantly, crucial 
information on quantum numbers is readily available. 

Despite such appealing features, due to its variational na- 
ture the method just described lacks built-in indicators of its 
own reliability. This drawback can be circumvented by re- 



11 



lying on unbiased techniques, such as ED in the basis, 
that are used in providing benchmarks to vaHdate the restric- 
tion to the NNVB manifold. As a first step in this direc- 
tion, in Fig. 1 1 we plot results for the relative difference be- 
tween the ground- state energy of the model Eq. (1) obtained 
by solving the GEP in the NNVB basis and from ED in the 
basis, (£;^NVB _ E^^)/E^^, for an AT = 24 site clus- 
ter. (Qualitatively similar results are obtained for a less sym- 
metric cluster with TV = 30 sites, not shown here. We note 
however that finite size effects in the NNVB energy per site 
seem to be surprisingly large, as shown in App. C. The qual- 
itative result regarding the region of best match with ED 
seems to be stable with system size however.) As expected, 
since long-range VBs are required in accounting for long- 
range spin correlations on two-dimensional lattices^^, E^^^^ 
compares poorly against Eq^ for couplings expected to sup- 
port magnetic phases (see sections I, II, and III). Conversely, 

(^NNVB _ £;ED)/^ED < 5%^ -^^ ^^^^^ 

diagram where magnetically-disordered phases are likely to 
be stabilized, cf. Fig. 2; in particular, for J3 = small rel- 
ative errors are observed for 0.2 < J 2 ^ 0.3, in agreement 
with Ref. 23. The fact that a small number of VB configura- 
tions {eight NNVB states, as opposed to 19873 in ED in the 
= {) subspace - lattice symmetries being exploited in both 
cases) is able to reproduce the GS energy in an extended re- 
gion of the parameter space up to a relative error that can be 
as small as ~ 1.5% constitutes good evidence that the NNVB 
subspace may be able to capture the low energy physics of the 
magnetically-disordered phases. Of course, to confirm this 
statement it is crucial to go beyond a simple energy -based cri- 
terion and to compare the nature of the correlations contained 
in this variational wave function and the exact one. In what 
follows we proceed to a thorough characterization of four-spin 
correlations. 




FIG. 12. (Color online) Four-spin correlations [Eq. (8)] on an = 
24 site cluster for ( J2, J3) = (0.5, 0.3), obtained from (a) ED in 
the basis and (b) by solving the GEP [Eq. (10)] in the NNVB 
subspace. The reference bond, in both panels, is indicated by the 
thick-black line. 



2. Four-Spin Correlations 

Real space correlations. We compute the four-spin con- 
nected correlation function Cijki as defined in Eq. (8) in the 
ED section, where z, j and k^l are pairs of NN sites (dimers) 
on the honeycomb lattice. Cijki is readily evaluated by ana- 
lyzing the loop structure in the transition graphs {Lpj\Lpi) for 
non-orthogonal NNVB states, in terms of which the lowest- 
energy state, solving the GEP, [Eq. (10)] is expressed (for 
technical details on how to compute expectation values for 
NNVB states see Ref. 63). 

We once more gauge the validity of our variational ap- 
proach and compare the so-obtained results for Cijki against 
those from ED for the A/" = 24 site cluster and {J2-,J?>) = 
(0.5, 0.3) in Fig. 12. Semi-quantitative agreement is observed, 
and interestingly the correlations obtained from the variational 
approach seem to be systematically smaller than those from 
ED, suggesting that the exclusion of longer-range VBs has 
the effect that VBC order is underestimated (see below). Re- 
garding the particular kind of VBC order that is stabilized, the 
two sets of data in Fig. 12 are consistent with both colum- 
nar (Read-Sachdev) and plaquette VBC order, although the 



particularly strong correlations at the shortest range (those in- 
volving the dimers closest to the reference bond) seemingly 
favor the later scenario (see the Appendix D), as vindicated in 
Refs. 9 and 23. We remark that evidence in favor of plaquette 
VBC order is also found from the histogram analysis in the 
framework of the effective QDIM presented in Sec. IV B. 

We take advantage of the substantially reduced dimension 
of the NNVB subspace and compute four-spin correlations for 
considerably larger clusters, comprising up to A" = 72 sites. 
The spatial dependence of Cijki is depicted in Fig. 13 for the 
N = 72 cluster and ( J2, J3) = (0.5, 0.3). An even stronger 
resemblance to the correlation pattern for "pure" columnar 
and plaquette states is observed than for the smaller cluster 
with A" = 24 sites [Fig. 12(b)], suggesting that the observed 
VBC pattern is not merely a finite size effect. This observa- 
tion is further corroborated by the data in Fig. 14(a), where 
we plot \Cijki\ as function of r (the distance between dimers 
i,j and k,l) also for A' = 72 and (J2, J3) = (0.5,0.3): 
\Cijki{r)\ decays slowly with r, consistent with saturation at 
large distances, and furthermore positive correlations are ap- 
proximately twice as large as negative ones, as expected for 
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pure columnar (Read-Sachdev) and plaquette VBC states. 




FIG. 13. (Color online) Four- spin correlations [Eq. (8)] on an = 
72 site cluster for (J2, J3) = (0.5,0.3), obtained by solving the 
GEP Eq. (10) in the NNVB subspace. The thickness of the bonds 
is proportional to Cijki and dark-blue (pale-red) indicates positive 
(negative) values. The reference bond is indicated by the thick-black 
line. 
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FIG. 14. (Color online) (a) \Cijki \ [Eq. (8)] as a function of distance 
r for the N = 72 site cluster and ( J2, J3) = (0.5, 0.3), obtained by 
solving the GEP Eq. (10) in the NNVB subspace; positive and nega- 
tive correlations are discriminated, (b) VBC structure factor per per 
number of bonds, ASpiaq./coi. /^b, as a function of the inverse system 
size N~-^ for (J2, J3) = (0.5, 0.3), as obtained from ED and from 
the variational approach in the NNVB basis. ^S'p^j^q /Coi. denotes the 
structure factor obtained by eliminating the strongest correlations en- 
circling the reference bond and linear fits (full lines; data for N — 18 
and = 36 are excluded from the fit) extrapolate the data to the 
thermodynamic limit. In both panels, dotted lines are only guides to 
the eye. 



Dimer structure factors. Reliable extrapolation to the ther- 
modynamic limit is achieved from the analysis of VBC struc- 



ture factors Soc (see Eq. 9). is expected to scale like 
A/N and thus the existence of a VBC phases is sig- 
naled by a finite value of the bond order parameter . In 
Fig. 14(b) we plot 5piaq./Coi./A^b and ^^i^q./coi./^b (A^b is 
the total number of bonds considered in the sum) as a func- 
tion of inverse system size, as obtained from both 5^ -ED and 
by solving Eq. (10), for (J2, J3) = (0.5,0.3). We first no- 
tice that, in agreement with our previous remark in connec- 
tion to Fig. 12, larger values for 5'piaq./Coi. and ^pi^q./coi. 
are obtained with S'^-ED and that the restriction to the NNVB 
manifold seemingly underestimates VBC order in the present 
case. Linear fits to the NNVB data shown in Fig. 14(b) 
(data for = 18 and for the less symmetric cluster with 
AT = 36 sites are discarded when fitting) yield the estimate 

^p?aq./Coi. ^ 0-1^' ^1^^^ ^^1^^ 9/^0 = expected 
for the pure plaquette VBC state (see the Appendix D). 



Finally, we analyze the strength of VBC order throughout 
the parameter space and in Fig. 15 we plot S^^^^ ^^^^^ /N\^ for 
the N = 72 site cluster as a function of J2, J3 G [0, 1]. Unlike 
what happens for smaller clusters, for which correlations mis- 
matching the sign structure of the columnar/plaquette VBC 
patterns are found in part of the parameter space, four-point 
correlations for the N = 72 cluster are always fully consistent 
with plaquette VBC order up to the point where, for given J3 
and increasing values of J2, one enters a regime (highlighted 
in Fig. 15) signaled by the occurrence of successive ground 
state level crossings (see Fig. 21 left panel for J3 = 0.3), 
likely associated with the breakdown of a description solely 
in terms of NNVB states (see related discussion in App. C). 
Maximal values of ^S'p^^^ ^^^^^ /N\^ are observed just before 
the first such level crossing occurs. For comparison we re- 
fer to the same quantity obtained by ED in the 5*^ using the 

= 24 cluster in the left panel of Fig. 9. For larger J3 values 
the agreement is quite remarkable, however as J3 is reduced 
to zero, the VBC correlations are significantly reduced in the 
ED approach compared to the NNVB results. 

On the other hand, the region signalled on Fig. 9 (right 
panel) by a strong staggered signal, does not appear to occur 
in a parameter range where the NNVB method can be safely 
used due to a rather modest variational energy (see Fig. 11) 
and the successive ground state level crossings (see Fig. 15 
and Fig. 21 left panel for J3 = 0.3). It may seem surpris- 
ing that a staggered VBC state fails to be naturally captured 
by the NNVB approach. However two distinct arguments can 
explain this paradoxical situation, (i) As mentionned in sec- 
tion IIIC, no definitive evidence supports the fact that the 
ground state is a singlet (non magnetic) state at the thermo- 
dynamic limit. In this respect a divergence of S'stag. is only 
a signal of the spatial symmetry breaking associated to it, 
but does not preclude the possibility of a (magnetic) nematic 
phase which would obviously be out of reach of the NNVB 
method, (ii) If the ground state is a dressed staggered (singlet) 
state, a good variational energy in the NNVB scheme may 
be hard to reach due to the structure of the staggered NNVB 
configuration (depicted in Fig. 4 (c)) : contrary to the case of 
the columnar or plaquette states (see Fig. 4 (a) and (b)), res- 
onating NNVB configurations needed to dress the pure VBC 
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FIG. 15. (Color online) S^y^^^q^y JN^ for N = 72, as obtained by 
solving the GEP Eq. (10) in the NNVB subspace, as a function of 
(-^2, Js)- The radius of the circles is proportional to ^S'p^^^q /Coi./^b- 
Numbers correspond to lO^S'pi^^q /cq^ /A^b- 



state and lower its variational energy involve long resonat- 
ing loops in the corresponding overlap diagram (see App. D) 
hence producing exponentially small overlaps and corrections 
to the bare VBC energy. In summary, this effective locking of 
the staggered NNVB configuration could make it difficult to 
emerge in the NNVB approach. A route to cure this issue and 
allow a more efficient relaxation to other VB configurations 
may be to include longer range dimer configurations assum- 
ing that the GEP remains numerically tractable. 



value problem. In practice the transformation is conveniently 
carried out by a truncated diagrammatic expansion, contain- 
ing only the most relevant terms. This derivation is pre- 
sented in App. B, and provides a rather simple quantum dimer 
model with single- and double-hexagon resonance and poten- 
tial terms, that reads 



eff ■ 




Interestingly, the parameters of the effective QDM (given in 
App. B) only depend on the ratio J2/(Ji + J3), so that the 
physics is invariant along simple isolines in the J2 — J3 plane 
within the QDM description. In retrospect - comparing to 
Fig. 15 — this invariance is also exhibited approximately by 
the NNVB approach, and within the boundaries of the VBC 
phase also to some extent in the 5^ -ED approach (cf. left 
panel of Fig. 9). 

The leading contribution of the QDM, with terms on a sin- 
gle hexagon, is a simple Hamiltonian of the Rokshar-Kivelson 
form^^, with off-diagonal and diagonal terms with amplitudes 
t and V respectively, that has been studied in great detail in 
Ref. 68. It turns out that the ratio V/t = 1/4 is fixed (in- 
dependent of J 2 and J3), for which case it was shown that 
the ground state is in the plaquette phase. Since our effective 
Hamiltonian reduces to this form when J2/( Ji + J3) = 3/8 
(i.e. the two-hexagon terms vanish), we expect that plaquette 
physics will occur in this region. However, it is not yet clear 
what will be the extent of this phase, since when we move 
away from this line, the 2-hexagon terms appearing in the ef- 
fective QDM might alter this behavior. 



1. Comparison with NNVB 



B. Exact Diagonalization of an Effective Quantum Dimer 
Model 

The NNVB approach used in the previous section requires 
an extensive numerical treatment of the non-orthogonality of 
the VB states. It is therefore numerically demanding and pre- 
cludes the use of efficient iterative algorithms, such as the 
Lanczos^'^ algorithm. In this respect, it would be desirable 
to base the study on a reliable (orthogonal) quantum dimer 
model (QDM) in order to significantly increase the accessible 
system sizes. 

Recently, a generic scheme for the derivation of QDMs 
from underlying Heisenberg Hamiltonians has been proposed 
in the context of two-dimensional frustrated antiferromag- 
nets^^'^^. This method aims to transform the generalized 
eigenvalue problem of the Heisenberg model in the short 
range valence bond basis (which was discussed in the pre- 
ceding subsection IV A 1) to an effective orthogonal eigen- 



In deriving the QDM in App. B we made several assump- 
tions and simplifications. First of all we neglected subdomi- 
nant terms in order to keep the model simple. This may cause 
some approximation errors that can be detected by comparing 
the QDM with NNVB approach. Secondly, the QDM is de- 
rived for an infinite lattice^^, which substantially improves the 
finite size scaling behavior. As shown in App. C the ground 
state energy converges to the same value as the NNVB ap- 
proach, however for the QDM the convergence is much faster, 
i.e. it has smaller finite size corrections. This validates the 
quantum dimer model, and also explains the finite size differ- 
ences between the QDM and the NNVB approach. 

A careful comparison between both approaches is pre- 
sented in Appendix C. In order to illustrate on a specific ex- 
ample how all these methods agree, we have computed the 
finite- size gap to the first singlet excitation with momentum 
K for (J2, J3) = (0.5,0.3). This set of parameters was 
chosen based on ED correlations computed for an A/" = 24 
cluster (with a grid spacing of 0.1), since it provides a strong 
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FIG. 16. (Color online) Finite-size scaling of the singlet gap between 
the ground- state and the lowest excitation with momentum K for 
J2 = 0.5 and J3 = 0.3, obtained by ED in the basis, the NNVB 
and the effective QDM approach. Note that the singlet gap matches 
quite well between the different techniques, despite the somewhat 
poor variational energy of the NNVB and QDM approach. 



tential VBC candidates (columnar versus plaquette) is real- 
ized. 




FIG. 17. (Color online) (a) all dimers belonging to the same of three 
possible columnar states have the same dimer vector associated. A 
resonating plaquette contains two different dimer vectors, that con- 
tribute equally to the resulting histogram, (b) the phase space built 
by the dimer vectors forms an equilateral triangle, where the corners 
represent the columnar states, while the plaquette states are signaled 
by a binomial distribution on the edges. 



plaquette structure factor ^vbc (see Eq. 9 for its definition). 
As has been discussed already, both plaquette and columnar 
VBC have the same discrete symmetry breaking, correspond- 
ing to a 3 -fold degeneracy in the thermodynamic limit. Eor in- 
creasingly larger cluster sizes, the lowest singlet excitation at 
the two equivalent K points should collapse onto the ground- 
state. Moreover, because of the finite correlation length in the 
VBC, this singlet gap must ultimately vanish exponentially 
with increasing system size. 

In Eig. 16, we plot the scaling of this singlet gap A com- 
puted with all our numerical techniques. Unbiased ED data 
already shows a clear indication of a vanishing gap in the ther- 
modynamic limit, but the scaling seems rather ~ 1/N since 
we cannot reach large enough cluster sizes. Using NNVB 
data, we can extend our computations to larger clusters (up 
io N — 96), and we observe an excellent agreement when 
comparing to ED data. This is not obvious since for instance 
the variational NNVB energy is not that accurate (see Eig. 1 1 
and App. C), but computing energy differences can give accu- 
rate results when there is a systematic deviation in all energies. 
Looking at the scaling of NNVB data, we see a behavior that 
could be compatible with a scaling faster than 1/A^, but there 
are still some irregularities in the finite- size effects due to dif- 
ferent cluster shapes. Our last set of data is obtained from the 
QDM model, that was simulated up to a A/" = 126 cluster: 
we note a semi-quantitative agreement with other techniques 
for the gap numerical data. Moreover, as explained above, the 
QDM has much weaker finite- size effects, which is clearly ob- 
served in the plot where finite- size scaling is much smoother. 
The possibility to access large clusters with small finite-size 
effects allows us to show convincingly that the gap collapses 
fast enough and that the system has long-range VBC order in 
the thermodynamic limit. 

However, it remains to be determined which of the two po- 



2. Dimer Vector Histograms 

Studying (orthogonal) quantum dimer models offers two 
advantages: first of all one can study larger systems. Here 
for example we were able to study honeycomb samples with 
up to = 126 sites, the second, even more interesting point 
is that one has access to new observables which are hard to 
define and implement in either the or the NNVB basis. In 
order to detect the underlying phase of a Hamiltonian, one 
usually measures correlations in the ground state, as done in 
Sec. IV A 2. However, the QDM allows for the computation 
of a related useful observable^^"^^ . The idea is to associate 
a two-dimensional vector to every dimer and collect a his- 
togram of the vector occupations. Writing the ground state 
IV^o) = as superposition of orthogonal dimer con- 

figurations one defines the appropriate histogram as 

where C^^^Ny is indexing all dimer states \ipi) that have a total 
dimer vector 

{N,,Ny)= V[,,,]. (13) 

The left panel of Eig. 17 illustrates a particular choice of 
dimer vectors, that assigns three different vectors to the three 
columnar (Read-Sachdev) states. The phase space of the 
resulting histogram forms a triangle illustrated in the right 
panel, where the corners of the triangle represent the colum- 
nar states, while the plaquette states are signaled by a bino- 
mial distribution on the edges of the triangle. A staggered 
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FIG. 18. (Color online) Normalized dimer histograms P{Nx, Ny), as defined in Eq.(12), obtained within the effective QDM for different 
system sizes at J3 = 0. We observe pronounced plaquette VBC signals at the larger J2 values shown, with a tendency towards a reduced 
radius and a more [/(I) -symmetric behavior as the Neel phase at smaller J2 is approached. 



VBC state on the other hand would contribute to the center of 
the histogram. 

In Fig. 18 we display dimer vector histograms obtained 
within the QDM approach for the three largest samples with 
A/" = 72, 96 and 126 sites, and for several J2 values at J3 = 0. 
The parameter region in which the QDM approach is ex- 
pected to be appropriate for the original spin model is high- 
lighted with a white background, whereas the remaining pa- 
rameter region is shaded in grey. Let us start the discussion 
at J2 = 0.4, which is close to the point J2 = 3/8 = 0.375, 
where the QDM reduces to the Rokhsar-Kivelson model at 
V/t = 1/4, expected to display a plaquette VBC ground 
state^^. Indeed the histogram displays pronounced peaks in 
the middle of the edges of the triangle, as expected for a pla- 
quette VBC phase. These results thus corroborate the comple- 
mentary findings based on four-point correlation functions us- 
ing ED in the basis and the NNVB approach. Pushing the 
QDM somewhat further into the "unphysical" region of larger 
J2, the plaquette signal is even more pronounced. On the 
other hand by lowering J 2 towards J2 ~ 0.2, the histogram 
becomes more rounded and fuzzy, reminiscent of the emer- 
gent U(l) symmetry at deconfined quantum critical points at 
Neel to VBC transitions^^'^^"^^ . Upon lowering the J 2 param- 
eter, the radius is also somewhat reduced and the dimer cor- 
relation length increases, however our effective QDM using 
nearest-neighbor valence bonds only is not able to reproduce 
the vanishing of the VBC order parameter (i.e. the radius of 
the distribution) as the Neel phase is approached. In spite of 
this limitation, the approximate /7(l)-like symmetry exhib- 



ited by the histograms upon approaching the Neel phase may 
well be a physical feature of the Neel to VBC transition on the 
honeycomb lattice. 

V. CONCLUSION 

In the present work we have analyzed the phase diagram 
of the frustrated Ji — J2 — J3 spin- 1/2 Heisenberg model 
on the honeycomb lattice by using a combination of differ- 
ent ED approaches and a SCMFT treatment. We have local- 
ized the boundaries of several magnetically ordered phases in 
the the region J2, J3 G [0, 1], and found a sizable magneti- 
cally disordered region in between. We characterize a large 
part of this magnetically disordered region as a plaquette va- 
lence bond crystal phase. Interestingly we find that a partic- 
ular parameter- free Gutz wilier projected tight-binding wave 
function has remarkably accurate energies compared to finite- 
size extrapolated ED energies along the transition line from 
the well-known Neel phase to the plaquette VBC, a fact that 
points to possibly interesting critical behavior - such as de- 
confined criticality - across the transition. In contrast a direct 
Neel to staggered VBC transition has recently been shown to 
be strongly first order^^. 

Compared to previous work on the Ji — J2 — J3 phase dia- 
gram we localize precisely the magnetic phases (phases I and 
II in the phase diagram shown in Fig. 2) which have been dis- 
cussed to be present at the semiclassical level^'^^'^^'^^"^^, and 
we discuss the possibility of a reentrant collinear magnetic 
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phase IV in a region at larger J2, which would nevertheless be 
compatible with the staggered dimer correlations found ear- 
lier in the relevant region^' The possibility of a plaquette 
phase has been discussed recently in Ref. 23 restricted to the 
J3 = line, while an earlier work^ reported that the dimer 
correlations might not be sufficiently strong for a plaquette 
VBC, and put forward the idea of an RVB liquid. 

Here we established that a plaquette phase does indeed 
occur for larger J3 values when leaving the Neel phase (I) 
by computing dimer-dimer correlations via exact diagonaliza- 
tions in the NNVB subspace and the analysis of dimer his- 
tograms within an effective QDM, and we showed that the 
phase has a sizable extent in the J3 direction, including the 
magnetically disordered region found recently on the J2 = J3 
line^^ (and thereby clarifying its nature). The precise fate of 
the plaquette VBC upon approaching the J3 = line is how- 
ever still an open question. 

The situation regarding the staggered VBC versus magnetic 
order in phases III and IV in Fig. 2 is not clear yet, and the 
possibility of incommensurate behavior of spin correlations or 
magnetic order renders an ED analysis quite challenging. It is 
likely that this question can be more meaningfully addressed 
using coupled-cluster or spin f-RG techniques, as done re- 
cently in the context of incommensurate spin correlations on 
the frustrated square lattice^^ . 

At the technical level it is notable that we have found an 
interesting example where we could explicitly show that it is 
possible to derive an effective quantum dimer model which ac- 
curately describes the magnetically disordered plaquette VBC 
region. Such a connection was conjectured to be present al- 
ready some time ago^^. However, no precise connection be- 
tween a QDM and original spin models could be made at 
that time. It is still an open question to understand why the 
NNVB and the QDM approach are currently unable to detect 
and describe the staggered VBC (lattice nematic) discussed 
previously. It might be that both methods are biased towards 
dealing with valence bond configurations which retain some 
flipability on short loops, while the staggered VBC configu- 
rations do not contain short flipable loops at all. On the other 
hand it could also be that the lattice nematic state actually is 
a magnetically ordered state (at least in some part of param- 
eter space) which breaks the same lattice symmetries as the 
the staggered VBC, giving rise to qualitatively similar dimer- 
dimer correlations. 

Returning to one of the initial motivations — the under- 
standing of the magnetism of the half filled Hubbard model 
upon lowering of U/t, and the possible explanation of the spin 
liquid behavior found in Ref. 15 — the questions are: i) what 
is the effect of the sub-leading next nearest neighbor J2 cor- 
rection to the nearest neighbor Ji Heisenberg interaction in 
terms of new phases arising; and ii) are the required values of 
J2 / Ji for new physics beyond the Neel phase reachable by 
downfolding the Hubbard model to a spin model at intermedi- 
ate U/t, or does one have to consider more correction terms? 

Regarding i): the scenario developed in the present paper is 
that J2 destabilizes Neel order somewhere between J2/J1 ~ 
0.17 — 0.22, and then a plaquette VBC phase (or a disordered 
version thereof) sets in, up to a value of J2/ Ji ~ 0.35 — 0.4. 



For even larger values of J2 / Ji a lattice nematic {staggered 
VBC) [c.f. Refs. 9, 18, and 23 and right panel of Fig. 9] as 
well as magnetically ordered spiral phases can arise. Based 
on the success of the Gutzwiller projected Dirac sea wave 
function to quantitatively describe the energies at the Neel 
to plaquette VBC transition, as well as the qualitative agree- 
ment with respect to spin- spin and dimer-dimer correlation 
functions, we suggest that a deconfined critical point scenario 
might describe this particular Neel to VBC transition on the 
honeycomb lattice. Our current ED tools are admittedly not 
perfectly suitable to resolve more complex scenarios, such as 
an SU(2) algebraic spin liquid region^^ with a small but finite 
extent. The same is true for a small Z2 spin liquid region^^"^^ 
appearing between the Neel and plaquette phase. We also note 
that a recent instanton analysis of one kind of Z2 spin liquid 
revealed an instability to a VBC phase^^, in agreement with 
the plaquette VBC phase we find. Our analysis is however 
at variance with the phase diagram put forward in the VMC 
study of Ref. 50. In that work the succession of phases is Neel, 
a rather large Z2 spin liquid region, followed by a rotational 
symmetry breaking state. 

Regarding the question ii) a recent estimate on the ratio of 
J2 / Ji in the spin liquid region of the honeycomb Hubbard 
model was put forward based on continuous unitary trans- 
formations in Ref. 84, and a value of about J2/J1 ~ 0.06 
was quoted. While this value seems to be almost sufficient 
to enter a new phase in the Ji — J2 model according to the 
VMC analysis of Ref. 50, which reported a critical value of 
J2/ Ji = 0.08, our ED based value for the critical ratio is at 
around three times as large (0.17 — 0.22). We currently be- 
lieve that the small J2 / Ji value in Ref. 50 is due to a compar- 
atively poor variational energy of the Neel state, when com- 
pared to our finite size extrapolated ED energies, therefore 
shifting the VMC transition to a too small J2 / Ji value. So we 
believe based on our results that a simple Ji — J2 spin model 
alone does not allow a quantitative description of the spin liq- 
uid phase discovered recently in the Hubbard model^^. More 
work is needed to understand whether the phase adjacent to 
Neel phase at J3 = is a plaquette VBC with a small order 
parameter or a genuine spin liquid phase, in which case the 
Ji — J2 model at small J2 would at least qualitatively explain 
the physics of the Hubbard model on the insulating side of the 
Mott transition. Future efforts will also have to explore the ef- 
fects of higher order corrections and thereby unravel whether 
a quantitative spin-only description of the spin liquid phase in 
the Hubbard model on the honeycomb lattice is possible. 

Note added: After submission of this work we became 
aware of Ref. 85, where a spin f-RG study of the same model 
is presented. In that paper a considerably large magnetically 
disordered phase is also found. 
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FIG. 19. (Color online) Gutz wilier projected Dirac sea: Finite size 
scaling behavior and extrapolation L ^ oo of the spin correlations 
relevant the energy of the frustrated honeycomb Hamiltonian (1). See 
text for the explanation of the two types of symbols shown. 



Appendix A: Properties of the Gutz wilier-Projected Wave 
Function 

Wave function approaches to strongly correlated systems 
can be valuable approximations because a good variational 
wave function can provide significant physical insight due to 
its relative simplicity. In the present study we focus on a (com- 
pletely) Gutz wilier projected half-filled nearest neighbor hop- 
ping tight binding model on the honeycomb lattice. As this 
corresponds to a filled Dirac sea we term the wave function 
Gutzwiller projected (GP) Dirac Sea^^. 

Here we use a standard Monte Carlo procedure to evaluate 
correlation functions of the parameter free wave function ac- 
cording to the update scheme proposed by Ceperley, Chester 
and Kalos^^. While doing so we noticed the occurrence of 
slowly equilibrating starting configurations which had a sig- 
nificant effect in determining some of the correlation func- 
tions and their error bars. It is presently not clear to us whether 
this is due to a inefficient Monte Carlo sampling or due to a 
fat tailed distribution for some observables, as for example 
discussed in Ref. 88 for continuum systems. 

First we determine the nearest neighbor, next-nearest and 
third-nearest neighbor spin-spin correlation function, as they 
allow us to determine the variational energy of this wave 
function for the J1 — J2 — J3 Heisenberg Hamiltonian [Eq. (1)] 



studied in this paper. The finite size expectation values for 
lattices with N = 2 x L'^ with L = 8, . . . , 14 are displayed 
in Fig. 19. Two different data sets are shown, first the bare 
estimates with error bars including all independent Markov 
chains (open symbols), and a second set where the anomalous 
samples were removed in calculating the mean and the error 
bars (hatched symbols). The later procedure yields estimates 
which show a markedly smoother finite size behavior and are 
used to linearly extrapolate the estimates to L ^ 00 The 
energy per site is then found to be approximately: 

E{Ji, J2, J3)/N ^ -0.353 X 3/2 x Ji 
+0.128 X 3 X J2 
-0.120 X 3/2 X J3 . (Al) 

As already shown in Fig. 6 (for J3 = in the left panel and 
J3 = 0.3 on the right panel), the energy per site of this wave 
function is very close to the extrapolated energy per site of 
the frustrated Heisenberg model close to the supposed Neel 
to plaquette transition. This success is quite striking for a 
parameter free wave function. 




FIG. 20. (Color online) Dimer correlations [Eq. A2] evaluated in 
the Gutzwiller projected Dirac sea for a L = 11 sample. The black 
bond denotes the reference bond, while the blue (red) bonds denote 
positive (negative) correlations. The width of the bonds are propor- 
tional to the value of the correlation function. The (E) letter indicates 
hexagons with all negative correlations, while (?) indicates plaque- 
ttes with a staggered signal. The four bonds indicated by arrows are 
the only ones that differ in sign from the expectations for a Read- 
Sachdev or plaquette VBC. 

This surprising observation raises the question whether the 
wave function exhibits appropriate correlations to describe 
such a transition beyond the competitive ground state energy. 
We have therefore determined the spin correlation functions at 
larger distances and found that correlations are perfectly stag- 
gered according to the Neel pattern and decay algebraically as 
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1/r^^ with a decay exponent (jg ~ 1.7(2). Next we have mea- 
sured dimer-dimer correlations functions of nearest neighbor 
bonds 



(A2) 



kinetic term and vq = tQa^ for the potential one. Here we 
choose the bipartite convention with a = 1/ >/2. 

Interestingly, the amplitudes and vq depend only on one 
parameter, Jf^ = J2 / {Ji -\- Js) • Note, that this qualitatively 
agrees with the phase diagrams (Figs. 2 and 3) suggested ear- 
lier in this paper. One can therefore simplify the Hamiltonian 
to 



J2 



and display the correlation pattern for an L = 11 system in 
Fig. 20. As already pointed out in Ref. 50, the short range 
structure of the dimer correlations in this wave function is sur- 
prisingly analog to the one found in the spin liquid phase of 
the honeycomb Hubbard model . What has however not been 
noticed previously is that beyond the four "inverted" bonds 
(highlighted by arrows in Fig. 20) all the other dimer corre- 
lations explored here match the signs expected for columnar 
(Read-Sachdev) or plaquette VBC states as derived in App. D. 
The dimer correlations also seem to decay as a power law, but 
we have not been able to determine the corresponding decay 
exponent accurately enough. 

The correlations measured in the Gutz wilier projected 
Dirac sea qualify this wave function as a viable candidate to 
describe a critical state separating a Neel ordered magnetic 
phase from a VBC of columnar (Read-Sachdev) or plaquette 
type. The fact that this wave function simultaneously exhibits 
staggered Neel fluctuations, as well as columnar/plaquette 
VBC fluctuations, is reminiscent of the SU(2) algebraic spin 
liquid state on the honeycomb lattice put forward by Her- 
mele^^, which is however believed to describe an extended 
spin liquid region, instead of a single critical point^^. Further 
work is required to understand whether this wave function 
could possibly also represent a deconfined quantum critical 

point^4,75 separating the two phases or whether there is indeed the fusion of two flips on a hexagon can be written as 



,0 . (B2) 



One can check easily, that this relation also holds for 
processes that connect dimer configurations defined on two 
hexagons. While those contributions cannot be obtained an- 
alytically, we find some iterative, numerical algorithm, that 
allows for calculating the amplitudes for all possible terms of 
this kind. This algorithm appears to converge rapidly and will 
be briefly described in the following. 

While it is relatively easy to obtain the inverse of an opera- 
tor within the present scheme^^, calculating the square root is 
much less obvious. We therefore have to go beyond previous 
works in order to derive an expression for The idea is 

to explicitly work in a basis that is formed by all the diagrams 
that are considered. Hence, it is possible to write every sum of 
processes as a vector and every fusion as a linear map applied 
to this vector. As an example, in the basis 




an extended algebraic spin liquid phase present between the 
two ordered phases (Neel - plaquette VBC) discussed in this 
paper^^. Yet a different scenario has recently been advocated 
in Refs. 50, 78-80, where a gapped Z2 spin liquid has been 
proposed as a phase neighboring the Neel ordered phase. Vari- 
ationally the Z2 spin liquid was found^^ to be a tiny fraction 
lower in energy than the GP Dirac sea studied above. Whether 
this is also true beyond the variational realm remains an open 
question. 



Appendix B: Derivation of an effective Quantum Dimer Model 

In this section we redefine the Ji-J2-J3-Heisenberg model 
on the honeycomb lattice (1) as 



1 \ /I 
10 0-0 
10 0/ \0 

resulting in a contribution for a potential term on a hexagon 
and a kinetic term on two hexagons. 

Generalizing and applying this procedure to a larger basis 
allows for an iterative solution of 



to obtain O . Putting the result into Bl, we arrive at the 
Hamiltonian Eq. (11), with coefficients given by 



(Bl) 



where H is the matrix introduced in section IV A 1, (9 is the 
overlap matrix for the NNVB basis states and N is the number 
of sites. 

On the honeycomb lattice there is only one elementary pro- 
cess that resonates between the two possible valence bond 
coverings on a hexagon. As shown in Ref. 66, this naturally 
leads to a potential term, counting the number of flippable pla- 
quettes. The exact amplitudes of both processes are shown to 
be given by U = -(6J2 - 3Ji - ZJz)ol^I(\ - ot^) for the 



1 hexagon 


2 hexagons 


[| (2^f-l)] 


[| (8^f-3)] 


-0.6 


tio -0.049218(5) 




^10 0.001562(9) 



Note that changes sign at J 2^ = 1/2, while tio and 
change sign at J 2^ = 3/8. The ratio of ^lo/^io does not de- 
pend on J although its analytical value is not known at the 
present stage. We note in passing that the model with only 
the most relevant tio term has been studied in the context of 



19 



supersolids of hardcore bosons on the triangular lattice^^'^^, 
whereas the model at tio = ^lo = corresponds to a particu- 
lar point of the Rokhsar-Kivelson model studied in Ref. 68. 

The QDM combines the advantages from both exact diago- 
nalizations in the basis which can be performed efficiently 
based on the Lanczos algorithm and from the NNVB approach 
which reduces the Hilbert space significantly through the re- 
striction to nearest neighbor VB states. This approach makes 
it possible to study honeycomb samples of up to 126 sites us- 
ing space group symmetries'^ . 

One drawback of both the NNVB and the effective quan- 
tum dimer model approach is that they do not presently allow 
to gauge the quality of the approximation with respect to the 
Heisenberg model within the methods themselves. One there- 
fore needs to compare energies or overlaps with exact diago- 
nalization data of the original Heisenberg model for smaller 
system sizes in order to locate the regions in the phase dia- 
gram where the NNVB approximation is valid. 



Appendix C: Comparison Between NNVB and QDM 
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FIG. 22. (Color online) Ground- state energies for various sizes vs J2 
(J3 = 0.3) obtained by different numerical techniques. Extrapola- 
tions to the thermodynamic limit (see text) are also plotted. Inset: for 
J3) = (0.2, 0.3), finite-size scaling of the ground-state energy 
obtained with ED of NNVB and QDM models. 
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FIG. 21. (Color online) Low-energy spectra for the frustrated model 
Eq. (1), for J3 = 0.3 and as a function of J2, obtained by: (a) solving 
the GEP [Eq. (10)] in the NNVB subspace and (b) performing EDs 
for the effective QDM derived in the Appendix B.^"^ In both panels, 
results have been obtained by diagonalization of an = 54 site 
cluster and energies are relative to the ground- state energy. Note that 
this plot only serves to compare the NNVB and the effective QDM 
approach on a technical level, because the J2 values considered here 
are beyond the domain of validity of these approaches for the original 
Heisenberg model. 



Although both the NNVB method discussed in Sec. IV A 1 
and the approach relying on EDs of an effective QDM (Ap- 
pendix B) are similar in spirit, for both of them are formulated 
in terms of NNVB degrees of freedom and are thus especially 
suitable to the study of quantum spin liquids and VBC states, 
they differ somewhat in detail. One such difference concerns 
the fact that, in deriving an effective QDM, the diagramatic 
expansion detailed in Appendix B must eventually be trun- 



cated, but one lacks built-in indicators of the convergence of 
the resulting expression. On the other hand, overlaps are ex- 
actly dealt with within the NNVB approach (Sec. IV A 1), 
which is therefore immune to this problem, but this advan- 
tage comes at the cost that the system sizes that can be ana- 
lyzed via the NNVB approach are more restricted than those 
that can be handled by diagonalizing effective QDMs. Fur- 
thermore, and as far as finite-size analysis is concerned, the 
formalism described in Appendix B has the advantage that 
the amplitudes appearing in the effective QDM are computed 
on an infinite lattice,^^ implying that faster convergence to the 
thermodynamic limit is attainable within this approach. Al- 
together, these features imply that the formalisms detailed in 
Sec. IV A 1 and Appendix B should be regarded as comple- 
mentary to one another. In this sense, extensive comparisons 
between the results obtained from both methods and, due to 
the variational nature of the NNVB subspace, from unbiased 
techniques such as ED are clearly called for. 

As a step toward this goal, in Fig. 21 we compare the 
low energy spectra obtained from NNVB and by diagonal- 
izing the effective QDM derived in the App. B for the spin 
model Eq. (1) with J3 = 0.3 and varying values of J2 (in 
both cases, an = 54 site cluster has been considered). We 
first remark that overall features are similar in both spectra, in 
spite of the subtlety that energy levels displaying similar de- 
pendence on J2 are characterized by different quantum num- 
bers in Fig. 21(a) and Fig. 21(b).^'^ Another feature salient in 
Fig. 21 concerns the fact that, by increasing the value of J2, 
one enters a regime characterized by the occurrence of succes- 
sive level crossings. Note that this plot only serves to compare 
the NNVB and the effective QDM approach on a technical 
level, because the J2 values considered here are beyond the 
domain of validity of these approaches for the original Heisen- 
berg model. 

We proceed to a more systematic comparison and in Fig. 22 
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we plot the ground state energy dependence on J2 for J3 = 
0.3, as obtained from NNVB and EDs for the effective QDM, 
for system sizes N = 24 and 42 (data for other N are shown 
only in the inset, but are fully consistent with the analysis 
that follows). We first notice that much stronger finite size 
effects are indeed observed for the NNVB data, in agreement 
with our discussion above. In extrapolating to the thermody- 
namic limit we heuristically assume that the scaling relation 
Eq/N ^ A^~2 , only justified in the case of the Neel phase 
(see Sec. IIIB), also applies in the present case. As shown 
in the inset in Fig. 22, this indeed seems to be the case. Ex- 
trapolated values for the ground- state energy computed from 
NNVB and from the analysis of the effective QDM are also 
plotted in Fig. 22, and from the excellent agreement obtained 
we conclude that the dominant terms are correctly taken into 
account by the truncated expansion detailed in Appendix B. 
Finally, extrapolated data from both approaches based upon 
NNVB states are compared against those from EDs in the 
basis: we observe that agreement is optimal around the region 
where plaquette VBC order is strongest for J3 = 0.3 (Fig. 15) 
and where a description based on NNVB states should be at 
its most accurate level. 



Appendix D: Correlations in pure VBC states 

In this Appendix, we compute the expectation values of 
the 4- spin correlation function for the four candidate VBC 
states denoted l^/^c) (Columnar), {ipst) (Staggered), \il)sw) (s- 
wave plaquette) and (d-wave plaquette) in the thermo- 

dynamic limit (see Fig. 4). For the plaquette state indeed, 
we may consider s-wave or d-wave linear combinations of the 
two VB coverings of a single hexagon. For infinite systems 
each of these states is degenerate since it breaks spatial sym- 
metries. This degeneracy is lifted at finite size and, in order to 
allow direct comparison with finite size numerical results, we 
consider symmetrized trial states with (0,0) momentum and 
belonging to the trivial point group representation Al. 

Orthogonality. The overlap between two distinct 

components of {ipa) vanishes exponentially. This point is 
rather obvious for lipc) and {i/jst) but deserves more attention 
for \^sw) and Generically (rj^i) = 2-^(^'^')-^/2 

with N the size of the system and ni{i^j) the number of 
loops of the overlap diagram obtained by superimposing the 
dimer coverings i and j. A direct inspection of such a dia- 
gram shows that ni{i^j) = N/6 for the columnar state and 
ni{i^j) = \/N/2 for the staggered state, hence showing that 
any pair of distinct components becomes orthogonal in the 
thermodynamic limit. 

The plaquette state cases a = sw and a = dw are slightly 
more involved since includes 2^/^ overlap contribu- 

tions (see Fig. 23). It is possible, albeit not very illuminating, 
to find an upper bound of this sum of terms that goes to zero 
when the systems size goes to infinity. 

In fact, such a result is general: the overlap between two 
periodic states {ip) and {ip') related by a discrete symmetry 
*S is either 1 or in the thermodynamic limit. Before actu- 
ally showing this result let us mention how it can be antici- 




FIG. 23. Overlap {ipl, IV^^) between two distinct plaquette state com- 
ponents IV^^) (blue) and (red). 



pated using a physical argument. The two states being peri- 
odic, the structure of the scalar product {^p\^p') is itself peri- 
odic. It is thus tempting to infer that ~ a^"" in the 
thermodynamic limit, where Nc is the number of local pat- 
terns (scaling like the number of sites) and a is related to a 
local overlap or fidelity. In this case, either |?/^) = and 
a = 1 or ^ and a < 1 which implies = 
for an infinite system. While qualitatively correct, the scal- 
ing ~ a^^ is actually non-trivial. Indeed, the scalar 
product {^p\^p') does not generically break into a product of 
local disconnected terms but may involve arbitrary scale reso- 
nances. 

Let us consider a tensor product state |?/^) = 0c I ^c). where 
the same structure | (^c) defined on a cluster c is repeated on the 
lattice. The state |?/;') is related to |?/^) by applying the unitary 
operator S. Typically in our case, c is a hexagon, and 5 is a 
translation that transforms a hexagon into a neighboring one. 
Denoting the density matrixp = \ip){ip\ = 0cPc , the overlap 
can be written. 



(V^IV^O = Tr = Tr (^Jpc) . (Dl) 

But since Spc and Sp'^ do not commute in general, the re- 
lation (V^IV^^) = He (^'^Pc^ ^o^s hold, which illustrates 
the point raised previously according to which cannot 
be interpreted as the product of local quantities. 

However using the Holder inequality for traces we have for 
any finite TV, 



1(^1^01 = |Tr (^c^Pc) 
< Tr mc^p, 



< 



n Tr 



Spc 



where |X| denotes (X^X)^/^ and Nc is the number of clus- 
ters c (scaling linearly with the system size N). 
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Taking to infinity, Nc also goes to infinity and 



where Xo{X) stands for the maximal eigenvalue of the 
positive- semidefinite operator X. It is then straightforward 
to obtain the inequahty, 



|(^|^')l< lim {fc\S\vc 

Nc^oo 



(D2) 



Two cases can occur : (i) \(pc) an eigenstate of S in which 



case \ 
impHes 



= 1 or (ii) \(pc) is not invariant under S which 

{^c\S\^c) <landKV^|V^OI=0. 



Correlations. Considering the bond permutation opera- 
tors P5, it is straightforward to remark that (?/^Jj,|P5|?/^^) and 
(t/^J^IP^P^/ 1?/^^) vanish exponentially to as well since these 
operators can only produce local reconfigurations of loops. 
It follows that the three components of lip a) generate in- 
dependent contributions to the 4-point correlation function 
{Pij Pki ) — {Pij)'^. Its expectation values for the four trial VBC 
states is depicted in Fig. 24. Note that {PijPki) — {Pij)'^ = 

4({(s,.s,.)(Sfe.so) - {s,.s,n 

Note that in Ref. 23, the authors claim that the three pla- 
quette states are not orthogonal in the thermodynamic limit, 
which is in contradiction with our general result. However, 
their approximate numerical values for the dimer-dimer corre- 
lations agree with our exact ones. In Ref. 9, the dimer-dimer 
correlation between parallel bonds on neighboring hexagons 
is quoted to be 0.01, while we find a negative value of —0.09, 
a result which agrees in sign with ED data deep in the plaque- 
tte phase (see Fig. 8c). 
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